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Single receptor-ligand bonds have finite lifetimes, so that biological systems can dynamically 
react to changes in their environment. In cell adhesion, adhesion bonds usually act cooperatively 
in adhesion clusters. Outside the cellular context, adhesion clusters can be probed quantitatively 
by attaching receptors and ligands to opposing surfaces. Here we present a detailed theoretical 
analysis of the stochastic dynamics of a cluster of parallel bonds under shared constant loading and 
with rebinding. Analytical solutions for the appropriate one-step master equation are presented for 
special cases, while the general case is treated with exact stochastic simulations. If the completely 
dissociated state is modeled as an absorbing boundary, mean cluster lifetime is finite and can be 
calculated exactly. We also present a detailed analysis of fluctuation effects and discuss various 
approximations to the full stochastic description. 

I. INTRODUCTION 

Cells in a multicellular organism adhere to each other and to the extracellular matrix through a large variety of 
different receptor-ligand bonds T]. Although not probed this way in traditional affinity experiments, adhesion bonds 
in physiological situations usually have to function under mechanical load. For example, cell-matrix adhesion in con- 
nective tissue is mainly provided by focal adhesions, which are based on transmembrane receptors from the integrin 
family connecting the actin cytoskeleton to the extracellular matrix. Focal adhesions of fibroblasts, the main cell 
type in connective tissue, are usually loaded by actomyosin contractility, in particular during tissue maintenance and 
wound healing. An important class of adhesion contacts in endothelial sheets are adherens junctions, which are based 
on transmembrane receptors from the cadherin family connecting the actin cytoskeletons of different cells. Endothe- 
lial tissue often is subjected to considerable external stress and strain, for example in lung and blood capillaries. 
Leukocytes circulating with the blood flow tether to and roll on vessel walls through transmembrane receptors from 
the selectin family connecting the actin cytoskeleton to carbohydrate ligands on the opposing surface. Here contact 
dissociation is accelerated due to the shear flow pulling on the cells. In general, there are many more physiological 
conditions in which adhesion clusters are subject to forces arising from intra- or extracellular processes, including cell 
motility, development and angiogenesis. 

During recent years, the behavior of different adhesion bonds under force has been investigated extensively on the 
level of single molecules by dynamic force spectroscopy |^ |^ |^ . This field has been pioneered by AFM-experiments 
by the Gaub group j5(] and later put onto a firm theoretical basis by Evans and Ritchie 6\. Because bond rupture 
can be modeled in the framework of Kramers theory as thermally assisted escape over one or several transition 
state barriers, bond strength is a dynamic quantity which depends on loading rate. Experimentally, this prediction 
has been impressively confirmed for different molecular systems 0i IM fej 10]- Dynamic force spectroscopy has been 
implemented with different experimental techniques, including atomic force microscopy 0, laser optical tweezers 
|8(| and the biomembrane force probe Hoj . The behavior of molecular bonds under force can also be probed in 
parallel plate flow chambers. Here usually the loading process is much faster than bond dissociation, which therefore 
effectively occurs under constant load 0, 0| . By now, dynamic force spectroscopy has shown that adhesion bonds 
feature a much more complicated behavior under force than suggested by the traditional affinity experiments in 
solution . Using concepts from the theory of stochastic dynamics 0, Q, 0, 0, 0| , a binding energy landscape 
can be reconstructed from the experimental data. During recent years, this has been accomplished for many different 
adhesion receptors, including integrins 18, 19], cadherins 2Qj and selectins Jl_, 22j. However, while dynamic force 
spectroscopy up to now has mainly been applied to single bonds, in physiological settings adhesion receptors usually 
operate cooperatively within clusters p^. Therefore the physical description of single adhesion bonds under force 
now has to be extended to clusters of adhesion bonds under force. Clusters also open up the possibility of rebinding of 
broken bonds, which is known to be essential to achieve physiological lifetimes of adhesion clusters. For single bonds, 
rebinding usually cannot be studied due to elastic recoil of the force transducer after bond rupture In contrast, 
for adhesion clusters open bonds can rebind as long as other bonds are closed, thus keeping the spatial proximity 
required for rebinding. Only if the completely dissociated state is reached, rebinding becomes impossible and the 
cluster disintegrates as a whole. 

Although it is clear that force leads to accelerated cluster dissociation, it is usually not known how it is distributed 
over the different closed bonds in different situations of interest. In many cases, most prominently in rolling adhesion, 
only few of the different bonds are loaded to an appreciable degree, thus dissociation occurs in a peeling fashion 
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|24l 125. 26]. However, due to geometrical reasons, even in this case there wih be a subset of bonds which are loaded 
to a similar extend. In the same vein, the loading situation at focal adhesions can also be expected to be rather 
complicated. For the case of homogeneous loading, one further has to distinguish between loading throu^ soft and 
stiff springs [23 ■ In the latter case, all bonds are equivalent and a mean field description can be applied 28]. In the 
first case, force is shared equally between all closed bonds and the coupling between the different bonds in the adhesion 
cluster is non-trivial. Recently, dynamic force spectroscopy has been applied to this case for the first time |29| |. Here, 
a vesicle functionalized with appropriate ligands is sucked into a micropipette and pressed onto a cell. On retraction, 
the vesicle is peeled off from the outside to the inside of the contact region. However, due to rotational symmetry 
around the micropipette axis, all bonds in a ring around the periphery of the contact area share the homogeneous 
loading. 

The equilibrium properties of adhesion clusters has been theoretically studied before [2M Isol IbH Is^ , mainly in 
reference to experiments on vesicle adhesion through specific ligand-receptor pairs |3^l34Ll35l |. For the non-equilibrium 
dissociation of adhesion clusters under force, a deterministic model has been introduced in a seminal paper by Bell 
[23. This model has been mainly used to study more specific problems, for example leukocyte rolling in shear flow 
Pd|. Recently, the deterministic Bell-model has also been extended to treat linear loading of a cluster of adhesion 
bonds, which usually is applied in dynamic force spectroscopy j23,l22|. A stochastic version of the Bell-model has 
been introduced, but studied only in the large system limit and for specific parameter valu es [ 371 . Later the stochastic 
model has been treated with reliability theory in the special case of vanishing rebinding [33|. Other special cases of 
the stochastic model have been treated in order to evaluate specific experiments, for example the binding probability 
between ligands and receptors on opposing surfaces as a function of contact time |3 US • 

In this paper, we use the stochastic version of the Bell-model to study the case of constant shared loading in 
comprehensive detail. In contrast to applications to specific experiments, we focus on generic features of the stochastic 
dynamics of a cluster of parallel bonds under shared constant loading and with rebinding. A short report on our main 
results has been given before As shown elsewhere, the same stochastic framework as used here for the case of 
constant loading can also be used to study the case of linear loading ji^l • Compared with the deterministic model, the 
stochastic model has several advantages: first, only the stochastic model allows to treat the experimental situation that 
rebinding becomes impossible once the completely dissociated state has been reached. Second, it includes fluctuations 
and non-linear effects, which are important for small adhesion clusters. Third, using the well-developed theory on 
master equations, the stochastic model allows to derive analytical results for cluster lifetime as a function of cluster 
size, rebinding rate and force, which are very helpful in evaluating adhesion experiments, including rolling adhesion 

In the following, we consider the situation in which a certain number of bound adhesion receptors has been clustered 
and connected to some force-bearing structure. We then ask how strongly the force accelerates dissociation, and in 
which sense dissociation can be balanced by rebinding. Since we are concerned with generic features of contact 
stability, our model does not consider spatial or concentration degrees of freedom. In Sec.m we define the stochastic 
variant of the deterministic Bell-model, which has three dimensionless parameters. Next we introduce the appropriate 
one-step master equation describing the stochastic dynamics of an adhesion cluster under shared constant force and 
with rebinding. We also explain how this master equation can be solved numerically with the Gillespie algorithm 
for exact stochastic simulations. In the two following sections, we discuss two special cases of the model in which 
considerable analytical progress can be made. In each case, we first discuss deterministic results, and then turn to 
the full stochastic model. In Sec. lIIII we discuss the case of vanishing rebinding. In this case, broken bonds cannot be 
reformed and the number of closed bonds in the adhesion clusters decreases in a unique sequence of rupture events. 
This can be used to construct a solution for the master equation and to derive an expression for the average lifetime 
of an adhesion clusters. In Sec. IIVI we discuss the case of vanishing force. In this case, we deal with a linear problem 
and analytical solutions of the master equation can be derived for a refiecting boundary. They can be used in turn 
to derive an approximation for the case with an absorbing boundary. Cluster lifetime can be calculated exactly as 
mean first passage time using Laplace techniques. In Sec. El we consider the general case with finite rebinding and 
finite force. Although full analytical solutions are only feasible in the case of small clusters, cluster lifetime can be 
calculated exactly for arbitrary cluster size. For larger clusters, full solutions of the master equation are obtained 
by exact stochastic simulations. Simulations are also essential to characterize single unbinding trajecories and to 
understand the role of fiuctuations. We close in Sec. IVII with a discussion of experimental issues. 
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FIG. 1: Schematic representation of an adhesion cluster under constant shared force: there are Nt = 5 receptor-hgand pairs, 
i = 3 of which are closed and equally share the load /. A closed bond ruptures with the dissociation rate k = koe-^^'. The 
Nt ~ i — 2 open bonds rebind with the force-independent eissociation rate kon- Our model has three parameters: cluster size 
Nt, dimensionless rebinding rate 7 = kon/ko and dimensionless force /. 



II. MASTER EQUATION 
A. Derivation 



The rupture of molecular bonds can be modelled in the framework of Kramers theory as thermally activated escape 
over a transition state barrier |^ 0, . Assuming an infinitely sharp transition state barrier leads to the so-called 
Bell equation for the single molecule dissociation rate as a function of force, k — k^e^/^'' [2^. Here the force scale 
Fb = ksT/xb is set by thermal energy kgT and the distance Xh between the potential minimum and the transition 
state barrier along the reaction coordinate of rupture. For a typical value Xb — 1 nm and physiological temperature 
T ~ 300 K, we find the typical force scale Fh — 4 pN. Physiological loading has indeed been found to be in the pN- 
range, both for cell-matrix adhesion |4^l43. l45| and rolling adhesion [lllli^ . Values for fco and Fb have been measured 
during recent years with dynamic force spectroscopy for different receptor-ligand systems, including integrins p5ll9j. 
cadherins [23| and selectins [21I . While the dissociation rate k depends mainly on the internal structure of a 
bond, the association rate kon includes the formation of an encouter complex and therefore depends on the details 
of the situation under consideration. It is very difficult to determine kon experimentally, especially in the case of 
cell adhesion, when the interacting molecules are anchored to opposing surfaces [s^. l39l 11^ . In order to focus on 
the generic features of cluster stability, here we assume that kon is a force-independent constant, in accordance with 
earlier theoretical work [23, 123, ISa, |33 ■ Future modelHng might refine this assumption, considering for example the 
effect of ligand-receptor separation controlled by polymeric tethers 143, "i^, "5^. 

For the following, it is convenient to use dimensionless quantities. We define dimensionless time r — k^t, dimen- 
sionless force / = F/Fb and dimensionless rebinding rate 7 = kon/ko. The dimensionless single molecule dissociation 
rate is fc/fco = e^. We consider a cluster with a constant number of Nt bonds, which initially are all closed and 
then undergo rupture and rebinding according to the appropriate rates. Since bond rupture is a discrete process, the 
stochastic dynamics of the bond cluster can be described by the one-step master equation jsij ] 

— = r,+ip,+i + gt^iPi-i - [ri + g,\pi , (1) 

where Piir) is the probability that i bonds are closed at time t. Here the and gi are the reverse and forward rates 
between the possible states i {0 < i < Nt). They follow from dissociation and association rates of single bonds as 



i"i — r{i) = ie-^/* and gi — g{i) — j{Nt — i) 



(2) 
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Our model has three parameters, namely cluster size iVj, rebinding rate 7 and force /. Since i > should be 
guaranteed at any time, tq — has to be set for / > 0, in addition to the definitions in Eq. Moreover, Eq. Q 
implies 50 > 0, that is, after rupture of the last closed bond new bonds are allowed to form. This corresponds to 
a reflecting boundary of the master equation at z = 0. As explained above, in biological and biomimetic situations 
rebinding of the completely dissociated state is usually prevented by elastic recoil of the transducer. Therefore in the 
following we set go = in order to model an absorbing boundary at j = 0. Because the values for ro and go do not 
follow the general form given in Eq. (|2Jl, the boundary at i = is an artificial boundary. Concerning the upper end 
of the set of states at i = Nt, the form gjVj = represents a reflecting boundary and guarantees i < Nf. Thus, the 
upper boundary is a natural boundary of the master equation. 

A quantity of large interest is the average number of closed bonds N{t) = (i) = X^iL'i *Pi (''')■ From the master 
equation Eq. one can derive ,52J 

f'E'S^ = -«■» + (*)). (3) 

If r{i) and g{i) were both linear functions in i, Eq. (PJ would become an ordinary differential equation for N. This 
suggests to study the deterministic equation 

^ - -r((z)) + 5((z)) = -iVe//^ + j{Nt - N) (4) 

as has been done by Bell |23|- Below we will see that the analysis of this equation gives valuable insight into the 
generic features of our model. However, it is important to note that for / > 0, the reverse rate r{i) in Eq. ^ is 
non-linear in i and the average in Eq. Q cannot be taken. Instead lower moments are related to higher moments and 
one arrives at a complicated hierarchy of coupled differential equations. The solution of the deterministic equation 
Eq. Q will therefore deviate from the average number of closed bonds obtained from the solution of the master 
equation Eq. Q. The same problem arises for the higher moments. For example, for the variance cr^ ~ {1^) — (i)^ 
one can derive [53 

^ = {9{^) + r{^)) + 2{{^ mgi^) m) , (5) 



where again the average cannot be taken. As an approximate treatment, one can expand r(i) in a Taylor series around 
the average for i, (i) = A'^ [5]| . Restricting the expansion to second order, thus assuming a Gaussian distribution, 
leads to the following equations 

^ - -NeJf^ + ^(N, -N)- -^e//^^ , (6) 



In principle, these equations can be solved by numerical integration. However, it is much more instructive to consider 
the original master equation. Moreover, the deterministic equation Eq. (0J and its improved version Eq. (|SJ) cannot 
describe the effect of an absorbing boundary at « = 0. In order to consider this experimentally relevant case, one has 
to study the master equation Eq. with the rates given in Eq. (0). Finally, only the full stochastic analysis reveals 
the detailed effect of fluctuations. 



B. Numerical solution 



Below we will present analytical solutions for several special cases of the master equation. In the general case, 
we numerically solve the master equation by Monte Carlo methods. In detail, for each set of parameter values Nt, 
f and 7, we generate between 10^ and 10^ trajectories with the help of the Gillespie algorithm for exact stochastic 
simulations |53. l54|. By averaging for given time r over the different simulation trajectories, we obtain the desired 
probability distributions {Pi(T)}jJg. In general it is also rather instructive to study single simulation trajectories, 
because their specific features are expected to be characteristic also for experimental trajectories. 

The Gillespie algorithm was originally developed for exact simulation of the stochastic dynamics of coupled chemical 
reactions. Applied to our case, open and closed bonds correspond to two different species of molecules and the 
transition between these two species, that is rupture and rebinding, correspond to chemical reactions. The Gillespie 
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algorithm is very efRcient because rather than discretizing time in small steps, it generates jumps between subsequent 
reactions. The basic quantity of the Gillespie algorithm is the probability P{ii,t\tq, X)dT that the next reaction 
occurs in the time interval [tq + r, tq + t + dr] and is of type fi under the condition that at time tq the system is in 
state X. In our case, fi has only two values corresponding to rupture and rebinding, and the state X of the system 
is completely described by the number of closed bonds i. Since the rupture and rebinding rates from Eq. Q are 
constant between subsequent events, P does not depend on absolute time tq. In fact it reads 

t|to, X) = Pi^i, t\i) = PoiT\i)a^ (8) 

where Pq is the probability that no reaction occurs in the time interval [0, r] and is the reaction rate for reaction 
fi. Pq satisfies the differential equation 

/ \ 

(9) 

and the initial condition Po{0) = 1, therefore Po{T\i) — e^(^''"^)^. P{ii,T\i) is properly normalised to unity as can 
be shown by integrating over time and summing over reactions. The Gillespie algorithm generates trajectories in 
which subsequent reactions are separated by the following rule. In the absence of other reactions, the probability for 
a reaction fi in the time interval [r, r + dr] is given by 

P^{t) = a^.e-'^-^dT . (10) 

The integral 

F,{r)^ r p,iT')dT' = 1-6-'^-^ (11) 








is the probability for a reaction occuring until time r. It increases strictly monotonically from to unity and thus 
can be inverted. In order to generate a random variable which is distributed according to Eq. H10|l . one generates 
a random number ^ which is uniformly distributed over the interval [0, 1] and inserts it into the formula 



HO 



(12) 



This is done for each type of reaction, leading to a set of times r^. The time for the next reaction is then chosen as the 
smallest r^, that is r = min^(Tp). As shown in Is^, '53|, this rule generates trajectories with the correct distribution 
of times and types of subsequent reactions. With the forward and reverse rates Eq. (0) for rebinding and unbinding 
of molecular bonds in the cluster, the random times are determined by the functions 

HO , HO , ^ 

^/ = - (T, and = . 13 

This algorithm is exact in the sense that the only sources of inaccuracy lie in the choice of the random number 
generator and the finite number of trajectories used to calculate probability distribution. 



III. VANISHING REBINDING 
A. Deterministic analysis 

We start our analysis with the case of vanishing rebinding, 7 = 0. Then the deterministic equation Eq. Q reads 

^ = _7Ve//^ . (14) 
ar 

In principle, the total number of bonds Nt is irrelevant in this case. However, it is reintroduced through the initial 
condition iV(0) = Nt. Then Eq. lO is solved imphcitly by (23| 

r(«)^E{i:\-E(r\ (15) 
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where E{z) = dz'e~^ j z' is the exponential integral Unfortunately, the inversion for Nij) is not possible in 
general. 

In the deterministic description, cluster lifetime T^et can be identified with the time r at which only one last bond 
exists. Setting N(Tdet) = f in Eq. ifTCjl gives 

T,,t^E(J^^ -E{f) . (16) 

From this result, we can extract three different scaling regimes. For small force, / < 1, we use the small argument 
expansion of the exponential integral, E{z) « — F — ln(z) (where F = 0.577 is the Euler constant), and find 

Tdet-\nNt. (17) 

This corresponds to the familiar case of radioactive decay, when the differential equation dN/dr ~ —N leads to 
exponential decay N — Nte~'^ . 

For intermediate force, 1 < / < iVt, we can rewrite Eq. 1)16(1 for the cluster lifetime as the sum of two integrals, 

e / , e 



Tdet = / dz + / dz . (18) 

Jf/N, z z 

The second integral can be estimated to be J-j^ dz{e~^ / z) < dze^^ — — e"^ < e^^. For the first integral, we 
can expand the integrand for small arguments, leading to /j/jVt dz{e^^/z) « If/Nt dz{l — z)/z « lii{Nt/f). Since 
Nt > f, the second integral can be neglected and we have 

Tdet « In (^y) . (19) 

For the time evolution of the cluster, both exponential integrals in Eq. H15|l can be replaced by the small argument 
approximation as long as N > f and hence the exponential decay proceeds until the force per bond is f/N = 1. 
Thereafter the decay will be faster than exponential due to the destabilizing effect of force. 

For large force, f > Nt, the second term in Eq. (|16l) can be neglected and we can use the large argument approxi- 
mation for the exponential integral, E{z) « e^^/(l + z), leading to f23j 

e-//JVt 

Therefore cluster lifetime decays faster than exponential with f/Nt in this regime. Now the small argument approxi- 
mation is not applicable to either of the two terms in Eq. (|15|1 and the decrease in N will be faster than exponential 
over the whole range of time. 

In summary, the analysis of the deterministic equation Eq. (|14|l allows to identify three scaling regimes of small, 
intermediate and large force. This analysis also shows that f/Nt is an important scaling variable, which we will 
therefore use to analyse also the stochastic case. 



B. Stochastic analysis 

For finite force, / > 0, the reverse rate r(i) in Eq. |(2Jl is non-linear in i and the boundary at i = is artificial. 
Therefore the master equation in general cannot be solved with standard techniques. However, in the case of vanishing 
rebinding, 7 = 0, one can use the fact that the decay of the cluster corresponds to a unique sequence of events, with 
the number of closed bonds decreasing monotonously from Nt to 0. The transition from state i (with i closed bonds 
present) to the state i — 1 (with one more broken bond) is a Poisson process with the time-independent rate . If i 
bonds are present at time t, the probability that the next bond ruptures at time t -I- r' is given by 

p,_,_i(r') = r,e-'-'^' . (21) 

The state probability pi-i is related to the state probability pi and the transition probability pi—>i-i by the recursive 
expression 

Pi-i{T) = / dr' pi{T')pi^i-i{T - t) ^ r^e^''''' dr' p^{T')e''^ , (22) 
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FIG. 2: The state probability Pi(r) is the probability that i bonds are closed at time t (0 < i < Nt). Here the pi axe plotted 
as a function of time r for a cluster of initial size Nt = 10 for 7 = and (a) / = 0, (b) / = 1, (c) / = 10 and (d) / = 50. 



which uses the fact that the state i — 1 can be reached only through the state i. This scheme is solved by 

Mr) = I n -(^■) E ^ ^'^'^^^ TT ZTTT^ \ (23) 
1 j=i+i / j=i 



Nt . 
r(j)T TT ^ 

Ji r(fc) - r(j) 



as can be shown by induction. The properties of this expression follow from the properties of the reverse rate 
r(i) = ie-^/^, which for finite force, / > 0, is a non-monotonous function of i. It diverges for i — > 0, has a minimum 
at i = [/] (the integer closest to /) and grows as i for z ^ oo. In the unlikely case that the value of / is such that 
r{j) = r(k) for j ^ k, the limit of the expression in Eq. (|23|) for r{k) r{j) has to be taken carefully. In order to 
treat this case properly, one has to replace Eq. (|23|l by 



Aft 



Nt 



Mr) = I n Kj) e 

ij=i+l I j=i 



Nt 

n 



( 



r{k) - r{j) 



J 



Nt 



n i 



k—i 
r{k)=r{j) I 

k^3 I ) 



(24) 



where we have used limA^o(l — e ^'^)/A = r for A = r{i) — r{k). In Fig.|21wc plot the full solution to the master 
equation, that is the state probabilities Pi(r) from Eq. if^ . for Nt = 10 and four different values of force /. For small 
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FIG. 3: Average number of closed bonds A*' as a function of time r for different cluster sizes Nt = 1,2,4,8 and 16. (a) 
f/Nt = 0.01, (b) f/Nt = 0.1, (c) f/Nt = 1.0 and (d) f/Nt = 5.0. 



force, all states are appreciably occupied during the decay, that is each of the curves is a maximum of the set of curves 
during a certain period of time. In the long run, po approaches unity and all other pi disappear, because without 
rebinding, the cluster has to dissociate eventually. For increasing force, the shape of the curves changes considerably. 
Now the lower states (with small number of closed bonds i) hardly become occupied during the decay process. For 
very large force, the maximum occupancy changes directly from the initial state Nt over to the detached state 0. 

With the help of the exact solution Eq. (|23J), any quantitiy of interest can now be calculated. One quantity of large 
interest is the average number N{t) = {i) — ^Pi^^) of closed bonds at time r. For a single bond, N is simply 

the probability pi that the bond is attached, 

7V(r)=pi(r) = e-^'^ . (25) 

For a two-bond cluster we have 

N{t) = pi(r) + 2p2(r) - {e^^- + (1 - eJ /^)e~^^"'^] . (26) 

For increasing Nt, the corresponding expressions become increasingly cumbersome. In general, N(t) is a sum of Nt 
exponentials with the different relaxation rates r{i) with 1 < i < Nt- For small force, / < 1, r{i) « i and the smallest 
rate corresponds to z = 1, that is, N ~ e""^ on large time scales. In this case, clusters of any size decay with the 
same slope as single bonds and the difference between single and multiple bond rupture lies in the prefactor, not in 
the time-scale of average decay. For intermediate force, 1 < f < Nt, and large time scales, decay is dominated by 
i = [/], that is ~ e~f^^ . Thus the absolute value of force governs the long time behavior, with different sizes 




FIG. 4: Average number of closed bonds A'^ as a function of time r for Nt = 8 and 16 and (a) f/Nt = 0.1 and (b) f /Nt — 1. 
Solid and dashed lines are stochastic and deterministic results, respectively. 




FIG. 5: Single simulation trajectories (solid lines) compared with the average number of closed bonds N{t) (dashed lines) for 
Nt — 10 and 25 and (a) f/Nt — 0.1 and (b) f/Nt — 1. In (b), in addition the deterministic results are plotted as dotted lines. 



showing up only in the prefactor. For large force, / > Nt, decay at large time scales is dominated by z = Nt, that 
is '--^ g-ATte^/^tr xhis implies that for a given force /, the largest clusters show the slowest decays in the long 
run. However, if one controls f/Nt rather than /, the cluster with the smallest size will decay the slowest, since it is 
subject to the smallest absolute force. In Fig. Owe plot logiV as a function of time r for different values of cluster 
size Nt and force per initial bond, f/Nt- All curves initially show an exponential decay with the rate of a single bond. 
For small forces decay stays exponential for almost all times. The larger force, the earlier decay crosses over to the 
late stage regime of super-exponential decay. 

As noted in Sec. |nl due to the non-linear form of r{i) = ie^ 1"^ for / > 0, the first moment N{t') of the stochastic 
solution Eq. (|23|1 is not identical with the function N(t) obtainted from the deterministic equation Eq. (|14|1 . In Fig. 01 
results for N(t) derived from the deterministic and the stochastic description are compared to each other. For a small 
but non-zero force, the non-linearity is small and the agreement between the two results is good in the initial phase of 
the decay. Towards the end of the decay strong deviations are observed. Here, the force on each bond grows strongly 
and the non-linearity of the transition rates is large. For increasing force, fluctuations become less relevant and the 
deviation between deterministic and stochastic results is increasingly restricted to the very end of the decay process. 

In Fig.|3jL and b the result for the mean number of closed bonds N(t) is compared to single simulation trajectories 
for small and large forces, respectively. The single simulation trajectories are expected to resemble experimental 
realizations for the time evolution of the number of closed bonds. The figure shows that for small force, the trajectories 
decay in a similar way as does the average. For large force, the trajectories decay in a more abrupt way than the 
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FIG. 6: Probability for dissociation of the whole cluster D as a function of time r for Nt — 1, 2, 4, 8 and 16 for (a) f /Nt — 0.0, 
(b) 0.1, (c) 1.0 and (d) 5.0. 



averages, that is they appear to run along the average for most of the time, but then decay rather abruptly towards 
the completely dissociated state. In this case, fluctuations do not so much affect the typical shape of the rupture 
trajectory, but rather the timepoint of rupture. The reason for this typical behavior is that a large fluctuation towards 
the absorbing boundary inevitably leads to a runaway process, since force is increasingly focused on less and less bonds 
due to shared loading. This type of rupture process is similar to avalanches or cascading failures in highly connected 
systems. Although rupture is rather abrupt, its timepoint is widely distributed, leading to the smooth decrease of 
N{t) observed in the average. In the large force case in Fig.jSJs, we also show the deterministic results for N{t) (for 
the small force case in Fig. they hardly differ from the stochastic results) . These curves show that the abrupt 
decay of single simulation trajectories at large force is somehow predicted by the deterministic description, compare 
Fig. 21 This had to be expected because the deterministic equation describes a representative yet single trajectory. 

The probability for dissociation of the overall cluster (that is for rupture of the last bond) is defined by D{t) — 
Po{t) = ripilr) and follows from Eq. (|23|) with the reverse rate ri from Eq. ||2Jl. The resulting formula has been given 
before in Ref. [s^. For a single bond it is simply D{t) = e^e~'^ . For Nt = 2 we have 

Dir) . (e-^^^ - e--^^^^) . (27) 

In the special case / = 21n2, the two rates r(l) and r(2) are equal and we have 

D{t) = IGre-^^ . (28) 

In general, as for N{t), D{t) is a sum of exponentials e"''*^*^'^ and the decrease on long time scales is governed by the 
exponential which decreases the slowest. In Fig.Elwe plot D{t) for different values of Nt and f /Nt (by controlling f /Nt 
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FIG. 7: Average adhesion cluster lifetime T as a function of f /Nt for cluster sizes Nt = 1, 10, 100, 1000 and 10000. The dashed 
horizontal curves are the harmonic numbers, which are good approximations in the small force regime, / < 1. The dotted 
curve is the approximation T = \n{Q.&lNt / f) for the intermediate force regime, 1 < f < Nt- 



rather than /, the curves have comparable averages). The case iVi = 1 is a Poisson process with simple exponential 
decay. For Nt > 1, D{0) = 0, because instantaneous rupture of all bonds at r = is a higher order process. For large 
times, all curves decay exponentially. For vanishing force, / = 0, the curves are very similar, with the same slope at 
large times. The maxima of the cluster dissociation rates for / = are described by T^ax = In A'^f, in agreement with 
the result Eq. H17|) from the deterministic description. For small f /Nt, the distributions are Gauss-like with small 
asymmetry and variance. For large f /Nt, they became Poisson-like, that is they develop a strong asymmetry with 
a maximum close to zero and a pronounced long-time tail. The reason is that in this case, decay is dominated by 
rupture of the first bond, that is we are effectively back to a single bond system (except that D{Q) as always for 
multiple bonds). 

The average cluster lifetime can in principle be calculated as the first moment of the overall dissociation rate 

/•OO 

T = / dr T D{t) . (29) 
Jo 

In practice, it has a simple form which can be derived without using the probability distribution Eq. H23f) . The waiting 
time spent in state i before the transition into state i — 1 is a stochastic variable characterised by the distribution 
function Eq. (|21|) . Its average is given by the inverse transition rate l/r{i). Since the decay process is a sequence of 
such independent Poisson processes, we simply have 

^ rr- (30) 
T = ^ - = (31) 

i=l 

which are the harmonic numbers. The lifetime of a two-bond cluster is increased by a factor 3/2 — 1.5 with respect 
to the single bond, that of the three-bond cluster by 11/6 — 1.8, and so on. For large Nt one can write [s^ 

r«lniVt + ^+r (32) 



E 



For / = we get |3a|5 
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FIG. 8: (a) Variance (j]v for the number of closed bonds in relation to the average number of closed bonds A'' as function of 
time r for f /Nt = 1.0 and Nt = 8 (solid) and 16 (dashed), (b) Variance ar for cluster lifetime in relation to average cluster 
lifetime T as a function of f /Nt for Nt = 1, 10, 100, 1000 and 10000. 



where T = 0.577 is the Euler constant. In fact this approximation is very good already for small values of Nt- The 
weak (logarithmic) dependence for large Nt means that for large adhesion clusters, size matters little since the bonds 
decay independently of each other and on the same time-scale. This result differs from the deterministic one for small 
force, Eq. H17(l . by the constant T and the additional contribution l/27Vt, which vanishes for large clusters. For small 
force, / < 1, Eq. (|32|l is a good approximation for cluster lifetime T. For intermediate force, 1 < f < Nt, the reverse 
rate grows rapidly for states with i < f, whereas for states with i > f, r(i) remains close to i. Therefore we can 
approximate the average lifetime of a cluster as Hj^^ — Hf. Using Eq. H32II . we get 

T^\n{Nt/f). (33) 

Thus cluster size Nt is now replaced by an effective size Nt/ f, as we have already found in the deterministic framework, 
compare Eq. H19|) . For large force, f > Nt, the only term which contributes to Eq. H3U|) is the one for the rupture of 
the first bond. Then 

(34) 

and we deal essentially with a single bond effect: if the first bond breaks, all remaining bonds follow within no time 
('domino effect'). This effect is also evident from the dissociation rate D{t), which for very large force approaches a 
Poisson distribution, compare Fig.|^. In Eq. H34|) . the numerator represents the probability for single bond rupture 
under force, while the denominator represents the probability that any one out of Nt identical bonds breaks first. 
Since f > Nt in this regime, the first effect dominates and T increases with Nt- For a given f /Nt, on the other hand, 
the lifetime decreases with increasing Nt, due to the increase in absolute force. In contrast to the deterministic result, 
Eq. H2U|) . the stochastic result Eq. (|34|) does not scale with f/Nt- In Fig. [T] we plot the average cluster lifetime T 
from Eq. H30|l as a function of f/Nt for different values of Nt- For small force, / < 1, T plateaus at the value given 
by the harmonic number iJjVj according to Eq. H31|l . In the regime of intermediate force, 1 < / < A^t, all curves fall 
on the master curve T = In (0.61(7V(//)), as predicted by Eq. H33|) . For large force, / > Nt, the scaling with f/Nt is 
lost, as predicted by Eq. H34|) . Although deterministic and stochastic predictions for cluster lifetime T have similar 
overall features, the deterministic result underestimates the plateau at small force and predicts an incorrect scaling 
with f/Nt at large force. 

Higher cumulants of the various distributions provide information about the effect of fluctuations. For the number 
of closed bonds at time r, the width of the distribution is described by the variance, defined by (J%{t) = {i^) — (i)^. 
In Fig. IH^i, we plot the relative standard deviation, aN{T)/N{T), for cluster sizes Nt = 8 and Nt = 16. It it zero 
initially due to the initial condition and diverges for large times. In regard to the distribution of cluster lifetime, 
the variance ax can be calculated in the same way as the average lifetime, because for a sequence of independent 
stochastic processes, all cumulants simply add up. The variance of the Poisson process Eq. (|21|l is l/r^(z). Therefore 
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the variance for cluster lifetime is 

^-T-E^y (35) 

For vanishing force this expression reads 

^T = E^=C(2)-^«(A^* + 1), (36) 
1=1 

where ip^^^{Nt + 1) is the trigamma function and the Riemannian ^-function. For increasing Nt, the variance 
converges to a finite value. For zero force this limit is given by 

4=E^=C(2) = ^, (37) 

i—l 

because the trigamma function vanishes in this limit. This result is an upper limit for the variance in general, because 
the reverse rate increases with increasing force, r(j) > i. The relative standard deviation ut/T of cluster dissociation 
is always smaller than unity, since (X)^)^ > X]^^- -'^'-'^ single bond rupture, we have a single Poisson process and 
it becomes exactly unity. For vanishing force and large clusters, it scales as ^ 1/lnA^t. Although it decreases with 

1/2 

increasing Nt, it does so in a different way than the Gauss process, which decreases as ~ 1/^t ■ The reason is 
that the contributions from the different subprocesses are not constant, but decrease as rupture proceeds. For large 
forces, f > Nt, cluster dissociation becomes a Poisson process governed by the rupture of the first bond. Then the 
first term dominates in Eq. (|30|l and Eq. H35|l . Therefore the relative standard deviation (Tt/T ~ 1 again. Moreover, 

1 /2 

now gt/T ^ l/-^t , because now only the first Nt — f subprocesses contribute, with roughly similar values, like in a 
Gauss-distribution. In Fig.|SjD, we plot (Tt/T as a function of f /Nt as it crosses over between the cases of vanishing and 
very large force, with a minimum around f /Nt ~ 0.3, that is in the intermediate force range. The narrow distribution 
at intermediate force is also evident in Fig. |H| Fig. also shows how the relative standard deviation decreases with 
increasing cluster size Nt. In general, the agreement between deterministic and stochastic descriptions is best for 
large cluster size Nt and intermediate force 1 < f < Nt- However, it should also be noted that the definition of 
deterministic lifetime is somehow arbitrary, because a discrete cutoff has to be introduced in a continuum description. 
Especially for small clusters the choice of the cluster size at which dissociation occurs will have a large influence on 
T. 



IV. VANISHING FORCE 

A. Deterministic analysis 

We now turn to the case of vanishing force, / = 0. Then the deterministic equation Eq. Q reads 

dN 

= -N + j{Nt - N) . (38) 

UT 

For the initial condition A^(0) = Nt, its solution is 

1-7)1" r 1 , , " 

Neq . (39) 



N{t) = li— Nt 

1 + 7 



1 + ie-(i+7)r 
7 



Thus there is an exponentially fast relaxation from N to the equilibrium state with Neq = jNt/{l + 7) closed bonds. 
Neq increases linearly with the rebinding constant 7 from Neq = for 7 = and saturates at Nt for 7 > 1. In the 
deterministic description, the lifetime of the cluster is infinite, because the completely dissociated state iV = is never 
reached. 



B. Stochastic analysis 



In the case / = 0, the reverse rates defined in Eq. ^ are linear in i and r(0) = at « = 0. Natural boundary 
conditions imply g(0) = 7A^t, that is a reflecting boundary condition at i = 0. A linear system with natural boundary 
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conditions can be solved with standard techniques. For the initial condition Pi{0) = Si^Nt^ a- generating function has 
been derived by Mc Quarrie [56| : 



A - /fs-lV-(i+'^)^ + l + '7s\^* 

G(.,r)=^.>.(r)^( (-- ^ (1+/^" ) • (40) 



=0 



The state probabilities follow from the generating function as 



1 d^G{s,T) 



(1 + 7)^* ■ ^ ^ 



One can easily check that by setting 7 = in Eq. (|41|l . one obtains the same result as by setting / = in Eq. H23|l . 
Eq. (|41|) shows that the systems relaxes to the stationary state on a dimensionless time scale 1/(1 + 7), thus the 
larger rebinding, the faster the system equilibrates. In the stationary state, the state probabilities follow a binomial 
distribution 

^(I')j^ (42) 

because the bonds are independent and each bond is closed and open with probabilities 7/(1 + 7) and 1/(1 + 7), 
respectively. 

The generating function also allows to calculate all moments of the distribution: 



9"G(s,r) 



a(lns)" 



(43) 



Since now r{i) is linear in i, the first moment N{t) = (i) is identical to the solution Eq. (|39|l of the deterministic 
equation. In order to assess the role of fluctuations, we calculate the variance: 



(l_e-(l+7)r) 
(1+7) 



^Ur) = {^') = ' , ' Nir) . (44) 



The relative standard deviation a^/N essentially scales as N for all times, thus fluctuation effects decrease with 
increasing bond number in the usual way. The stationary state value is lim^^oo aN{T)/N{T) = ((1 + 7)A^eg)^^/^ = 
(7A't)~^/^. Therefore larger rebinding does not only increase the equilibrium number of bonds, but also decreases the 
size of the fluctuations around N^q- This leads to a narrow distribution for large cluster under strong rebinding, with 
a small probability of coming close to the lower boundary. 

In Fig.Et., we plot the state probabilities pi from Eq. (|41(l for cluster size Nt = 10 and rebinding constant 7 = 1. 
The system quickly relaxes to the equilibrium state. The only difference for different initial conditions is in the initial 
transient. In particular, for A^o — N^q, the average does not change in time, although the distribution initially spreads 
to the binomial one. For 7 = 1, the stationary distribution is symmetric around the average. The width of the 
distribution for different 7 is illustrated in Fig. ^jp, which shows the average number of closed bonds normalised by 
the equilibrium number of bonds, that is N/N^q, together with the relative standard deviation, apf /N, for different 
values of the rebinding constant 7. The curves for iV are independent of Nt due to the normalization. Fig. |5Jd shows 
that with increasing 7, relaxation becomes faster and the width of the distribution decreases. 

For the biologically important case of an absorbing boundary at z = 0, it seems to be rather difficult to find a 
closed-form analytical solution for arbitrary cluster sizes. For the case Nt = 2, we will present such a solution in 
the next section. For arbitrary Nt, we use Monte Carlo simulations as described in Sec. |nj In Fig. ^1 we show 
individual simulation trajectories for different parameter values of interest, in comparision to the average number of 
closed bonds for reflecting and absorbing boundaries at z = 0. The plots show that the number of closed bonds in 
a cluster first relaxes towards the steady state value, for which rupture and rebinding balance each other. Although 
for the absorbing boundary the number of closed bonds decreases with time in average, for individual realizations it 
stays roughly constant, until a large fluctuation towards the absorbing boundary leads to loss of this realization. The 
time-scale for the decrease in N is thus determined by the probability for fluctuations from the steady state to the 
absorbing boundary. 

Because a full analytical solution is not available for the case of an absorbing boundary, we now introduce an 
approximation for this case. It is similar to the local thermal equilibrium description introduced by Zwanzig for 
modelling protein folding dynamics jSTj]. Our starting point is that for large clusters and strong rebinding, the 
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FIG. 9: (a) State probabilities Pi(r) from Eq. II41II for Nt = 10 with / = and 7 = 1 for a reflecting boundary at i = 0. (b) 
N/Neg and aN/N for Nt = 100, / = and 7 = 0.1, 1, 5 and 10. 




10 20 30 40 ^ 50 



FIG. 10: Single simulation trajectories for / = 0, 7 = 0.5, Nt = 10, 100 and 1000 and an absorbing boundary at i = 0. Solid 
lines are the average number of closed bonds A*" and dashed lines are the equilibrium number of closed bonds N^q. 



absorbing boundary is a small perturbation to the solution for the reflecting boundary, Eq. (|41|l . which in the following 
we will denote by {pi}^g. Since 5(0) = for the absorbing boundary, po = ripi with n = I and probability will 
only accumulate in the completely dissociated state. Since po is slaved to the other state probabilities and since we 
expect only a small perturbation for the states with i > 1, we assume that here the different boundary only leads to 
a simple renormalization caused by the 'leakage' into the absorbing boundary: 

K(r) =p,(r)(l-po(r)) for i>l (45) 

Po(r)= rpi{T')dT' . 
Jo 
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Since relaxation to the steady state is faster than decay to the absorbing boundary, Pii^r) can be taken to be the 
stationary value, that is the constant pi(oo) according to Eq. (|42|l . Then po{t) — pi(oo) Jq — po{T'))dT' , which is 
solved by Po{t) = 1 — e~P^^°°^^ . Therefore Eq. (|^ simplifies to 

Pi(T) =pi(oo)e-Pi(°°)^ for i>l (46) 
Po(t) = l-e-w(oo)r _ 

We conclude that the solution decays exponentially on the time scale l/pi(oo). In Fig. ^2 we plot Monte Carlo 
solutions for the state probabilities in comparison to the approximation. For j — 1, the approximation does not 
work well for Nt = 5, but it does so already for Nt = 10. For 7 = 5, the approximation works well for both cluster 
sizes. Note that in this approximation, a term po(oo)e~^^^°°''^ is missing for proper normalization ^fjloPi = 1- This 
is a small error for large clusters and strong rebinding. In order to assess the validity of Eq. (|46|l . we note that it 
presupposes that the time scale for relaxation to the steady state, 1/(1 + 7), is smaller than the time scale for decay 
to the absorbing boundary, l/pi(c>o) = (1 + j)^* /jNt- Therefore 7 should be larger than {NeqY^^* — 1. 

It follows from Eq. I|46|l that the mean number of closed bonds decay in an exponential way, N{t) — iVeqe"^^^""^'^. 
This is confirmed by Fig. 112b .. which shows the corresponding simulation results. For Nt — 2,5 and 10, we have 
Pi(oo) = 0.5,0.16 and 9.7 x 10~^. Numerically we find 0.6,0.13 and 0.01, thus the approximation is rather good. 
In Fig. 112b . we plot numerical results for the standard deviation (Tat. The initial increase of is well described 
by Eq. H44|) for the reflecting boundary, thus the boundary has little influence here. Large clusters stay close to the 
steady state during the time shown and the approximation is applicable. For small clusters, the variance grows larger 
than the steady state value before is decreases exponentially while the cluster size N approaches zero. The variance 
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FIG. 12: (a) Average number of closed bonds A'^ obtained from stochastic simulations of the master equation for 7=1 and 
Nt — 1, 2, 5, 10 and 100. (b) Variance crjv of the cluster size distribution for the same parameters as in (a). 



contains two time-scales. The second moment of the distribution decreases on the same timescale as the average, while 
the square of the first moment decreases twice as fast. The long time exponential decrease of (Tn is thus described by 
twice the relaxation time as is was found for the average number of bonds. 

Although an exact solution for the state probabilities seems to be impossible for the case of an absorbing boundary, 
more analytical progress can be made if one is only interested in the probability that the cluster dissociates as a 
whole. For the absorbing boundary, the cluster dissociation rate has been denoted by D{t) before. For the reflecting 
boundary, D{t) can be identified as the probability that the state i = is reached for the first time at time r if the 
system has started in the state i = Nt at time t — 0. This is a first passage problem which can be treated with Laplace 
techniques. Since the transition rates do not depend on absolute time, one can decompose the state probability for 
i = into two parts: 

P0(t)= r D{T')pofi{T^T')dT' (47) 



^0 

where po.o(''") is the state probability for state i = with initial condition Pi{Q) — Sifi. Po.q(t) can also be interpreted 
as the probability for having returned to the boundary after time t. A Laplace transform of the equation leads to an 
algebraic relation between the Laplace transforms of the three functions: 

Dis) = . (48) 

Here D{s) — e^'^'^ D{t) denotes the Laplace transform of the function D{t). The explicit form of the probability 
Po[t) is given in Eq. l|lT|l : 

P.{r) = ( j . (49) 

The probability po,o {t) can also be calculated with standard techniques [5^ : 

Po,o(r)= — ^ . (50) 



1 + 7 

The Laplace transforms of the these two functions are given by 



(1 + 7) ^ V « / 5 + HI + 7) 
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FIG. 13: (a) Average lifetime T of adhesion clusters as function of cluster size Nt for rebinding rates 7 = 0.0, 0.1, 0.5, 1, 5 and 
10.0. (b) Average lifetime T as function of rebinding constant 7 for cluster size Nt — 1, 2, 5, 10, 20 and 50. 



and 



Po.oGs)^ ^^ ^.7' , , (52) 

(1 + 7) V « / s + «(l + 7) 

so that the Laplace transformed first passage probability time distribution is 

Y^Nt fNt\ (-1)' 
2^i=0 V i ) s+i(i+7) 

Unfortunately, the analytical backtransform for D{s) seems to be impossible. However, the mean first passage time 
can be extracted from this result, because it does not require the backtransfrom [s^ : 



J. ^ dD[s) 



ds 



N, 



This equation is a polynomial of order A'^t — 1 in 7. The zero order term is the harmonic number H^t, so for 7 = 
we recover the result from Eq. (|31|l . In Fig. El we plot Eq. (|54ll as a function of cluster size Nt and rebinding rate 
7. As long as 7 < 1, cluster lifetime grows only weakly (logarithmically) with cluster size (at least for not too large 
clusters). For 7 > 1, the higher order terms in 7 take over and the increase in T becomes effectively exponential, as 
shown in Fig. 113b . In Fig. 113b . it is shown explicitly that increasing 7 to values larger than unity leads to a strong 
increase in lifetime. This effect is larger for larger clusters since the number of rebinding events in the dissociation 
path is larger. 



V. FINITE FORCE AND FINITE REBINDING 



A. Deterministic analysis 



Force destabilizes the cluster, while rebinding stabilizes it again. It has been shown by Bell that in the framework 
of the deterministic equation Eq. Q), the cluster remains stable up to a critical force fc |2j|. For the following it is 
helpful to revisit his stability analysis. In equilibrium we have 

iVe,e^/^- = 7(7Vt - A^eg) . (55) 

At small force /, this equation has two roots, with the larger one corresponding to a stable equilibrium. As force 
increases, a saddle-node bifurcation occurs. Above the critical force, no roots exist and the cluster becomes unstable. 
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FIG. 14: (a) Critical force Eq. 1571 in relation to total cluster size, fc/Nt, as function of the rebinding constant 7. It scales 
linearly (dotted curve) at small and logarithmically (dashed curve) for larger 7. (b) Stable steady value for number of closed 
bonds normalised by the critical value, Neq/Neq,c, for 7 — 0.1, 1 and 10 as function of force ///c for Nt = 100. Numerical 
results (solid lines) are compared to the approximation Eq. 1581 1 (dashed lines). The horizontal line marks the smallest possible 
value at the critical force fc. 



Exactly at critical loading, the two roots collapse and the slopes of the two terms become equal. This gives an 
additional equation 

e^^/^=(l-^) = -7. (56) 
These two equations allow to determine the critical values for cluster size and force: 

fc = Nt plog (^) and Nc = Nt ^ , (57) 

Ve/ l + plog(^) 

where the product logarithm plog(a) is defined as the solution x of xe^ = a. For small forces, the unstable fixed point 
is very close to zero. This implies that the stable fixed point is an attractor for most initial conditions. Close to the 
critical force, the unstable fixed point is close to the stable one and only the initial conditions above Nc will reach 
the stable fixed point. Eq. H57|l scales in a trivial way with Nt, but in a complicated way with 7. For 7 < 1, we have 
fc ~ ^Nt/e. Thus the critical force vanishes with 7, because the cluster decays by itself with no rebinding. For 7 > 1 
and up to 7 « 100, we have fc ~ 0.5Nt In 7. This weak dependence on 7 shows that the single bond force scale set by 
Ff, also determines the force scale on which the cluster as a whole disintegrates. Fig. 114b shows how fc crosses over 
from linear to logarithmic scaling with 7. 

Eq. (|55|l is an implicit equation which cannot be inverted to give Ncq as a function of the model parameters Nt, f 
and 7. In general, Ncq decreases from 7iVt/(l + 7) for / = to A^c for fc- For small forces we can find an approximate 
solution by first expanding the exponential function in Eq. I|55f) to second order in f/N^q and then expanding the 
resulting quadratic function for Ncq to second order in f /jNt: 



Ncq 



jNt 
1 + 7 



( f \ 7+1 / / 
\lNtJ 2 \jNt 



(58) 



Fig. 114b shows numerical results for Ncq/Nc in comparison with the low force approximation Eq. I|58|l for different 
rebinding constants 7 = 0.1, 1.0 and 10 as function of force f / fc and for cluster size Nt — 100. 

In the deterministic framework, cluster lifetime is infinite for f < fc, because a stationary state exists at Ncq. For 
/ > fc, cluster lifetime is finite, but strongly varies as a function oi Nt, f and 7. For / 3> fc, we can neglect rebinding 
and use the results from Sec. IIIII where we found for cluster lifetime 
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FIG. 15: Deterministic mean lifetime T of a cluster of Nt = 100 bonds for 7 = 1 as a function of the force-size ratio f/Nt- 
Numerical integration of the deterministic equation (solid line) is compared to the exact solution for vanishing rebinding (dashed 
line) and the inverse linear scaling (dotted line) predicted close to the critical force. 



compare Eq. H20|l . As force / is decreased from above towards the critical value fc, rebinding becomes important 
again and cluster lifetime diverges. To understand this limit, we note that here the system will evolve very slowly, 
because it is still close to a steady state. Therefore we can expand the time derivative of TV, compare Eq. Q), for 
small deviations from the critical state: 

In this limit, the lifetime will be dominated by the time spent close to the critical state. The time for a significant 
change A7V ~ — 1 in iV is 

T « Ar « e-(P'°s('^/°)+^' — - — . (61) 

f ^ fc 

Therefore T diverges like the inverse of f — fc- Fig. shows the lifetime of an adhesion cluster derived from 
numerical integration of the deterministic equation for 7 = 1 and Nt = 100 as a function of the force-size ratio f /Nf. 
The numerical results are compared to the approximation Eq. (|59|) for large forces and Eq. (|61|l for the divergence 
close to the critical point. Obviously both approximations work well for their respective limits. For different cluster 
sizes Nt, the plot remains basically unchanged (not shown), because the forces above fc are already in the range 
where Eq. 116|l predicts scaling with f/Nt alone. For different rebinding constants 7 the results are qualitatively the 
same, only that the critical force is shifted to different values. 



dN ^ d (dN\ 



B. Stochastic analysis 

In general, it seems to be difficult to find a closed- form analytical solution for the state probabilities Piir) for 
general values of 7, / and Nt. In the following, we will derive such an analytical solution for the case Nt = 2 with an 
absorbing boundary. In principle, solutions can be constructed in the same way for a reflecting boundary or larger 
clusters, but for increasing cluster size, the analytical procedure quickly becomes intractable. For this reason, we will 
later use simulations to deal with the general case. 
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We start by rewriting the master equation Eq. in matrix form: 

p — W ■ p . 

For the case — 2 with an absorbing boundary, p = {pq,Pi,P2)'^ a-nd 

/O ri 
W=i0 -in+gi) r2 
\0 gi -Ti, 



Eq. inil) is solved by |52 



p(T) = e^-.p(0) = ^c;,e^^. 



PA 



(62) 



(63) 



(64) 



where A and p\ are the eigenvalues and eigenvectors of W ^ respectively. The coefficients c\ have to be determined 
from the initial condition 



(o)-E 



c\P\ 



(65) 



Since the absorbing state po = (1)0,0)'^ is a stationary state, the corresponding eigenvalue Ag = 0. The other two 
eigenvalues are negative and correspond to transient states: 



Ai,2 = ~{^±uj) 



with r2 and lo being defined as 



f2 = (ri + r2 + 5i)/2 and lo = \/^^ — r\ri 
Note that < a; < and hence Ai_2 < 0. The transient eigenstates are 



1 



^A2 + ri 



Pi = — Ai + r2 and p2 



1 

51 



A2 + r2 
31 



(66) 
(67) 

(68) 



The three eigenstates p\ are linearly independent and form a basis of the state space of the cluster. With the initial 
condition pi = Si^fft i that is p(0) = (0, 0, 1)'^, the coefficients c\ follow from Eq. as 



A2 + r2 
A2 — Ai 



The final result then can be written as 

Po(t) 



= 1 - 



r2 



and c\. 



Ai + r2 
A2 — Ai 



(69) 



cosh(cjr) H sinh(wr 



Pi{t) — — sinh(a;r)e 



(70) 



cosh(cjr) 



smh(ti)r) 



There is a competition between the hyperbolic terms, which grow on the time-scale 1/w, and the exponential terms, 
which decrease on the time-scale l/f2. Since uj < fl, the exponential terms will win and only the stationary state 
survives. 

With the exact solution Eq. (|70|l . one now can calculate any quantity of interest. For example, the mean number 
of bonds, N{t) = J2i Wi^ follows as 



N{t) 



2 cosh(ajT) 



smh(ciJT) 



The dissocation rate for the cluster as a whole as given by D[t) — ripi — po, resulting in 
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FIG. 16: State probabilities pi as a function of time r for Nt — 10 with 7=1 for force / — 0.1, 1, 3 and 10 (a - d). In (a), the 
numerical solution is compared with the leakage approximation (dotted lines). The intermediate force values are chosen below 
and above the critical force fc = 0.278Nt. 



One easily checks that normalization is correct, D{T)dT = 1. Mean cluster lifetime T now follows as 

T^j TD{T)dT = -[2e'f + e-^'^ +-ie-^f/'^j . (73) 

As shown in the preceding sections for special cases, force leads to exponentially decreased lifetimes, while rebinding 
leads to polynomial terms in 7 up to order Nt — 1. 

Although the eigenvalue analysis can be used also for the general case of arbitrary cluster size, in this case it is 
more efficient to use exact stochastic simulations as described in Sec. Ill Bl In Fig. ^| numerical solutions of the state 
occupancy probabilities {pi}^o ^'^^ plotted for Nt = 10 with 7=1 for four different forces / = 0.1, 1, 3 and 50. This 
figure corresponds to Fig. [5] for vanishing rebinding and Fig. 1111 for vanishing force. For small force, the numerical 
solutions compare well with the approximation Eq. (|46|l introduced for vanishing force. For larger force, but still below 
the critical force, the state probabilities still decrease exponentially for large times, but the approximation Eq. H46|l 
breaks down, because the reference distribution {pi(oo)}jJo "^^^ had to be replaced by the unknown steady state for 
the case of finite force. If force is increased beyond the critical force {fc = 2.78 for 7 = 1), a simple description is 
not available, because equilibration and decay occur on the same timescale. For very large force, the behavior of the 
adhesion cluster approaches that for vanishing rebinding, with the analytical solution Eq. H23|l . 

Fig. 1171 demonstrates that the decay process changes dramatically as force is increased above the critical value. It 
displays trajectories of individual clusters with initially Nt = Nq = 10, 100 and 1000 closed bonds in comparison with 
the average number of bonds derived from a large number of these trajectories. Since fc — Q.278Nt for 7 = 1, Fig. 117b 
with / = 0.257Vt is below the critical value. For the largest cluster, the system equilibrates towards the steady state 
and then fluctuate around this value with very rare encounters of the absorbing boundary. For the smaller clusters. 
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FIG. 17; Single simulation trajectories for Nt — 10, 100 and 100 for 7=1 and at two different forces (a) / = 0.25Nt < fc and 
(b) f — Nt > fc- Representative trajectories are compared to the average number A'^ of closed bonds resulting from averaging 
over a large number of such trajectories. 
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FIG. 18: Average number of closed bonds N for Nt = 1, 10, 100 and 1000, 7 = 1 and (a) f/Nt = 0.25 and (b) f/Nt = 1. In 
(a), the numerical results are compared to exponentially decaying curves ~ e~'"^ (dashed lines) with a = 1.28,0.52,0.072 and 
0.0009 for increasing bond number. 



however, fluctuations towards the absorbing boundary frequently lead to loss of individiual realizations. As a result, 
the average number of closed bonds decays exponentially on a much faster timescale. Fig. II7b with f = Nt is above 
the critical force and the behavior is changed qualitatively. A steady state does not exist anymore and the clusters 
do not decay by fluctuations, but the size of each adhesion cluster is continuously reduced. Clusters of different size 
now decay on the same timescale and rebinding events are very rare in comparison to rupture events. 

Fig. ^1 plots numerical results for the average number of closed bonds N as function of time t for two different 
values of f/Nt and for cluster sizes Nt — 1, 10, 10^ and 10'^. For / — 0.25Nt < fc, after initial relaxation all curves 
decay exponentially. For f — Nt > fc, the larger clusters show a steep decrease in average cluster size at the end of 
the decay due to the effects of shared loading. For the small clusters, the average cluster size decreases slowly since 
cooperative effects are small. 

Fig. El plots the variance un(t) of the distribution pi for the two force values used in the two previous figures. 
Below the critical force the behavior is similar to that for vanishing force depicted in Fig. E] The variance decreases 
exponentially after having traversed a maximum. For forces above the critical force, a different behavior arises. After 
growing as expected in the initial phase, the variance displays a sharp peak. This effects becomes more pronounced 
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FIG. 19: Variance ajv for forces (a) / = 0.25iVt and (b) / = A^t for 7 = 1 and Nt = 1, 10, 100 and 1000. 





FIG. 20: Comparison of stochsatic and deterministic results for the average number of closed bonds A'^ derived from numerical 
solutions of the master equation (solid lines) and from integration of the deterministic equation (dashed curves) for cluster sizes 
Nt = 10, 100 and 1000 and forces (a) f/Nt = 0.1 and (b) 1.0. The rebinding rate is 7 = 1. 



the larger cluster size. 

In Fig. |2D1 a comparison of the average number of closed bonds in the stochastic and the determinsitic description 
is shown. N{t) is plotted for cluster sizes Nt = 10, 100 and 1000 for the forces f/Nt = 0.1 and 1.0 and the rebinding 
rate 7 = 1.0. For small forces / < /c, the average number of closed bonds equilibrates towards the steady state 
and remains constant thereafter. The fluctuations occuring in the stochastic description lead to a slow decrease of 
N. Above the critical force, the deterministic clusters decay as well, and in a more abrupt way than the stochastic 
average. 

We now turn to the dissociation rate of the overall cluster as a function of the model parameters. For 7 = 1 and 
/ = 0.25A^f and Nf, that is below and above the critical force, numerical results are plotted in Fig. [2] For a single 
bond, dissociation is a Poisson process with the maximum at r = and an exponentially decreasing dissociation 
rate D — r(l)pi = e-'e^'^ . For larger clusters and below the critical force, fluctuations to the absorbing boundary 
determine the rate of dissociation, which vanishes at r = 0, goes through a maximum and then decreases exponentially 
with time. As explained above, the exponential decay follows because decay proceeds by rare fluctuations from the 
steady state towards the absorbing boundary. Above the critical force, the dissociation rate for Nt > 1 becomes more 
sharply peaked and cannot be described with single exponential curves. A steady state does not exist anymore and 
dissociation does not proceed by fluctuations. The trajectories in Fig. El have shown that adhesion clusters decay 
fairly abrupt towards the end of the decay as a consequence of shared loading. This cooperative instability is the 
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FIG. 21; Dissociation rate D of the overall cluster for 7 = 1 and A^t = 1, 5, 10, 25, 100 and 1000. (a) / = 0.25Nt and (b) / = iVt. 



reason for the sharp dissociation distribution for large clusters under super-critical loading. The single bond that 
lacks these cooperativity, still shows the exponential dissociation rate which is now the slowest decaying for the given 
force size ratio. 

Whereas results for the dissociation rate have to be obtained numerically, the average lifetime can be calculated 
analytically [5ll |. The basic idea here is to sum the average times for any possible pathway leading from the initial 
cluster size A'^o towards dissociation at the absorbing boundary i = with its appropriate statistical weight. One can 
show that the lifetime Tf^^^No of a cluster with a total of Nt molecular bonds of which Nq are closed initially satisfies 
the equation [sj 

g{No) - TNt,No) + r{No) {Tn^^No~i ~ Tn^,n„) = -1 • (74) 

The left hand side can be considered to be the adjoint operator of the master equation acting on the average lifetime 
TNt,No- Foi' the initial condition A^o = ^t, the equation is solved exactly by 

Nt -, Nt~l Nt T-ri-1 />N 

T = Tj,,^,^y4T+y y il^^i^i^, (75) 

i=l ^ ' i=\ j=i+l llk=j-i ' 

where the first term is the result Eq. (|30|l for vanishing rebinding and the second term results in a polynomial of 
order A'^t — 1 in 7. For / = 0, Eq. 175|l is identical to the earlier result Eq. (|54|l obtained by Laplace transforms. Both 
expressions are polynomials of order A^t — 1 in 7, but in the general case from Eq. 175|) . the coefficients depend on 
force. For Nt — 2, we obtain the result from Eq. 173|) . For A^t = 3, we find 

r = e-^ + ^- + ^+7(^— ^+e-3//2j +72__ . (76) 

For Nt — 2 and 3, T can also be derived by explicitly summing over all possible dissociation paths. For larger Nt, 
direct summation becomes intractable and the results following from the general formula Eq. 175|) become rather 
lengthy. In general, force always affects most strongly those terms of highest order in 7, thus for 7 > 1, application 
of force is therefore an efficient way to reduce average lifetime T. For 7 < 1, T is dominated by those terms of lowest 
order in 7, thus here the reduction of lifetime with increasing force is not modulated by rebinding. 

Fig. 1221 shows the average lifetime of adhesion clusters of size A'o — 1, 2, 5, 10, 15 and 25 as a function of force-size 
ratio f /Nt for the rebinding constants 7 = 0.1 and 7 = 1.0. For small forces, / < 1, the average lifetime plateaus at 
the value given by Eq. (|54|l . For large forces, / > Nt, that is, when the force on each single bond is larger than the 
intrinsic force scale, the limit of vanishing rebinding applies (for Nt = 1 lifetime is independent of 7, compare also 
Fig. I15|l . The critical forces for the given rebinding rates are fc — 0.0355A^t and fc = Q.278Nt. In the intermediate 
force range, roughly around fc, the lifetime is reduced from the zero force to the zero rebinding limit. This reduction 
is dramatic for large clusters {Nt > 10) with appreciable rebinding (7 > 1), where the lifetime is reduced by orders 
of magnitude. We also show the lifetime following from the deterministic framework, which provides a lower limit 
for the lifetime at large forces, because here the largest clusters have the shortest lifetimes for a given force size ratio 
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FIG. 22: Average lifetime T according to Eq. 17511 (solid lines) of adhesion clusters with Nt = 1, 2, 5, 10, 15 and 25 as a function 
of ,f /Nt for (a) 7 = 0.1 and (b) 7 = 1. The critical forces for these rebinding rates are fc/Nt = 0.0355 and 0.278, respectively, 
where the deterministic results for the lifetimes (dashed lines) diverge. 




FIG. 23: Average cluster lifetime T (a) as function of rebinding rate 7 for Nt = 10 and f /Nt = 0, 0.1, 0.3, 0.6 and 1 and (b) as 
function of cluster size Nt for 7 = 1.0 and / = 0, 1, 2 and 10. 



f /Nt. Below the critical force the deterministic lifetime is infinite and the stochastic curves approach the plateaus 
Eq. 154|) determined by fluctuations towards the absorbing boundary. 

Fig. \2'6h demonstrates the influence of rebinding on the average lifetime at different levels of force. Here we show 
average lifetime T as function of 7 for Nt — 10 and for increasing values of force. For / = the curves are as depicted 
in Fig. 1131 Increasing force reduces the lifetime strongly and leads to an almost constant lifetime for different 7 
(compared the strong increase for / = 0). Only when rebinding is sufficiently strong that force is smaller than the 
critical force, / < /c, lifetime begins to grow. The increase observed then is similar to that for vanishing force, only 
that the absolute value of lifetime is smaller. For example, for / = 0.6Nt, the cluster grows strongly for g > 5 where 
the critical force is fc — Q.82Nt; for f = Nt the strong increase is observed for g > 10, for which the critical force is 
fc = l.lSA^t. A similar effect is observed for the dependence of average lifetime on cluster size, see Fig. 123b. At small 
Nt , cluster lifetime grows strongly at large forces according to Eq. (|30|) due to shared loading. For larger Nt , lifetime 
grows slowly until Nt is large enough that /c > / is reached. Above this size, T grows on a rate comparable to that 
for vanishing or small force. For 7 = 0.1, the increase of T with Nt is slow throughout the shown range of Nt- 
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VI. DISCUSSION 

In this paper, we have presented a detailed analysis of the stochastic dynamics of an adhesion cluster of size Nt 
under shared loading / and with rebinding rate 7. The corresponding master equation has been solved exactly for 
several special cases. For vanishing rebinding (7 = 0), the exact solution Eq. (|24ll could be constructed because cluster 
decay is a sequence of Poisson processes. For vanishing force (/ = 0), we deal with a linear problem, which can be 
treated with standard techniques. In the case of natural boundaries (that is for a reflecting boundary at i = 0), the 
exact solution Eq. 141|l follows with the help of a generating function. In the general case of finite force / and finite 
rebinding rate 7, for the case Nt — 2 and an absorbing boundary we used an eigenvalue analysis to derived the exact 
solution Eq. (|70|l . In principle, the same method can also be applied for a reflecting boundary or for larger clusters, 
but this does not lead to simple analytical results. 

For vanishing force (/ — 0) and an absorbing boundary at i = 0, we introduced the 'leakage approximation' 
(also known as 'local thermal equilibrium description' in the theory of protein folding), which treats the absorbing 
boundary as a small perturbation to the exactly solved case of the refiecting boundary. The resulting formulae given 
in Eq. 1)46(1 work well if average cluster lifetime T is much larger than the internal time scale 1/(1 + 7) (that is for 
large clusters or strong rebinding). All other cases have been treated with exact stochastic simulations using the 
Gillespie algorithm, which for large clusters is more efficient than the eigenvalue analysis. Moreover, the study of 
single simulation trajectories offers valuable insight into the typical nature of unbinding trajectories expected for 
experiments. 

Once the master equation is solved, either exactly or numerically, all quantities of interest can be calculated. In 
this paper, we focused on the mean number of closed bonds as a function of time, N{t), and the dissociation rate for 
the overall cluster, D{t). The first moment of D{t) then gives the mean cluster lifetime T . In this paper, we derived 
an exact solution T = T{Nt, f,'^) from the adjoint master equation, see Eq. H75|) . For the special cases of vanishing 
rebinding and vanishing force, we also showed how the exact formulae for T can be derived via completely different 
routes. The result for T = T{Nt, f) from Eq. H3U|) follows from the unique dissociation path without rebinding, while 
the result for T = T{Nt, 7) from Eq. (|54|l can be derived with Laplace techniques as a mean first passage time for the 
case of a reflecting boundary. In order to assess the role of fluctuations, we also calculated the standard deviations 
ctat and ctt for the distributions of the number of closed bonds and cluster lifetimes, respectively. 

A special focus of this paper was a detailed comparision between the stochastic and determinstic treatment. Re- 
garding mean cluster lifetime, the deterministic treatment is rather good in the case of vanishing rebinding, although 
it underestimates the plateau value for cluster lifetime at small force. In the presence of rebinding, the deterministic 
treatment fails, because it includes neither the effect of fluctuations nor the effect of an absorbing boundary. In 
particular, the deterministic treatment does not predict finite lifetime below the critical force /c, when clusters decay 
due to fluctuations towards the absorbing boundary. Only at very large force, when rebinding becomes irrelevant, 
does the deterministic treatment work well again. Regarding the average number of closed bonds, the deterministic 
model fails because it does not correctly treat the non-linearity in the rupture rate. This effect is most evident for 
small clusters and at late stage of rupture. In general, the mean number of closed bonds in the stochastic model decay 
in a smoother way than in the deterministic model, which typically shows an abrupt decay in late stage. This abrupt 
decay in fact is typical for shared loading and shows up in the stochastic model when one studies single simulation 
trajectories. In this sense, the deterministic model makes an interesting prediction which should be conflrmed in 
experiment, albeit not on the level of the first moment, as suggested by the deterministic model, but rather on the 
level of single trajectories, as suggested by the stochastic model. 

Our results can now be used to evaluate a large range of different experimental situations. The stochastic dynamics 
of adhesion clusters under force can be quantitatively studied with many different techniques, including atomic force 
microscopy, optical tweezers, magnetic tweezers, the biomembrane force probe, flow chambers, and the surface force 
apparatus. In all of these cases, by measuring cluster lifetime T and two out of the three parameters Nt, f and 7, 
the third parameter can be estimated with the help of our exact results. In general, our exact results nicely show 
how mean cluster lifetime T varies with cluster size Nt, force / and rebinding rate 7. For example, if the single bond 
lifetime was one second (fco = 1/s); for / = and 7 = a cluster lifetime T of one minute could only be achieved with 
10^^ bonds, because in this case, cluster lifetime scales only logarithmically with cluster size. However, for a rebinding 
rate 7 = 1 {kon = ^o): only Nt = 10 bonds are necessary, because lifetime scales strongly with rebinding, T ^ 7^'^^. 
Increasing force to / = 10 would decrease lifetime to T = 0.05 s, because T is exponentially decreased by /. To 
reach one minute again, cluster size or rebinding rate had to be increased such that f < fc- This implies Nt > 50 or 
7 > 10. It is important to note that these predictions are based on the assumption of rigid force transducers. In many 
experimental situations of interest, the force transducer will be subject to elastic deformations or even to viscous 
relaxation processes, like for example when pulling on cells |59( . In order to focus on generic aspects of adhesion 
clusters, here we only studied the minimal model for stochastic dynamics under force. 

Our results can also be applied to experiments in cell adhesion. For example, the biomembrane force probe 
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with linear loading has recently been used to study the decay of a^/Js-integrin clusters induced on the surface of 
endothelial cells K 

one makes sure that the clusters do not actively grow during the time of dissociation, 
similar experiments could now be done also for constant loading. Because in these kinds of experiments the exact 
cluster size is usually unknown, one had to convolute our results with a Poisson distribution for an estimated average 
number of bonds 34„ .39] ■ Recently, our result for the average cluster lifetime of two bonds under shared force and with 
rebinding, Eq. (|73(l . has been applied to the analysis of flow chamber data on leukocyte tethering through L-selectin 
[425. Since in this case force can be calculated as a function of shear flow, our formula can be used to the estimate 
rebinding rate, which in this case turns out to be surprisingly large. This in turn explains why dissociation dynamics 
in L-selectin mediated leukocyte tethering appears to be first order: for large rebinding, the leakage approximation is 
rather good, and decay is exponential. 

Our results can not be directly applied to adhesion clusters which compensate for force-induced decay by active 
growth, as it has been found experimentally for focal adhesions [g^I- Yet there are also interesting lessons for focal 
adhesions which can be learned from our model. For example, our stochastic analysis confirms the prediction from 
the deterministic stability analysis that cluster stability changes strongly around the critial value fc (although small 
clusters tend to decay also at smaller force due to fluctuations towards the absorbing boundary). It is interesting to 
note that recent experiments measuring internally generated force at single focal adhesions suggest that f/Nt, the 
most important scaling variable of our analysis, is roughly constant for different cell types |43l l45|. It is therefore 
tempting to speculate that focal adhesions (or subsets of focal adhesions) are regulated to be loaded close to the 
critical value fc/Nt — plog(7/e) from Eq. H57|l . In this way, cells could quickly increase force on single bonds by small 
changes in actomyosin contractility. Large force on single closed bonds in turn might trigger certain signaling events 
in focal adhesions, possibly by mechanically opening up certain signaling domains |6lj . Our speculation provides a 
simple way to estimate the rebinding rate, which is very hard to measure experimentally. Using compliant substrates, 
it has been found that focal adhesions are characterized by a stress constant ^ 5.5 nN//iTO^ 14^1^- We do not 
know which of the many different proteins in focal adhesions defines the weak link which most likely ruptures under 
force, but we expect that it will have a similar area density as the integrin receptors, which are expected to have a 
typical distance between 10 and 30 nm, corresponding to 10* and 10^ molecules per /im^, respectively. To obtain a 
lower estimate for 7, we therefore use Fc — 5.5 uN and Nt = 10*. For activated as/^i-integrin binding to fibronectin, 
recent single molecule experiments obtained for the molecular parameter values ko — 0.012 Hz and Fi, — 9 pN [l^ . 
Therefore the rebinding rate can be estimated to be at least 7 = 0.2, that is fco„ = 0.002 Hz in dimensional units. 
Based on future experimental input, it would be interesting to extend our model of passive decay to active processes 
resulting in cluster growth under force. 

Finally we want to comment that our model might also be applied to situations in materials science which are 
not directly related to biomolecular receptor-ligand pairs. One example is sliding friction, which recently has been 
modeled as dynamic formation and rupture of bonds under force [6^ . In general, we expect that many more cohesion 
phenomena in materials can be successfully modeled as dynamic interplay between rupture and rebinding. 
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